###  Author: David Huh
##
###  Dependencies: lmtest, sandwich
##
## generic function for calculating clustered sandwich estimates
sndwch <- function(model,family,rank,id,vcv,estfn,small.n) {
  
  ## degree of freedom adjustment robust variance estimator
  dfc.adj <- function(M,N,K,dist,smalln) {
    dfc0 <- (M/(M-1)) * ((N-1)/(N-K))  # gaussian (Stata)
    dfc1 <- M/(M-1)                    # non-gaussian (Stata)
    dfc2 <- M/(M-K)                    # Mancl & DeRouen
    
    ## abort if the small sample degree of correction is negative
    if (smalln & dfc2 <=0)
      stop("Small sample correction unavailable since number of parameters exceeds clusters")
    
    ## return df correction
    if (smalln) return(dfc2) else
      if (dist=="gaussian") return(dfc0) else return(dfc1)
  }
  
  N.obs <- length(id)            # number of observations
  N.clust <- length(unique(id))  # number of clusters
  
  dfc <- dfc.adj(M=N.clust, N=N.obs, K=rank, dist=family, smalln=small.n)
  
  if (missing(vcv)) vcv <- vcov(model)
  if (missing(estfn)) estfn <- estfun(model)
  
  uj <- apply(estfn, 2, function(x) tapply(x, id, sum))
  
  bread <- N.obs * vcv
  meat  <- crossprod(uj)/N.obs
  
  vcovCL <- dfc*(1/N.obs * (bread %*% meat %*% bread))
  
  return(vcovCL)
}